In [1]:
%matplotlib inline
#%matplotlib widget
from astropy.cosmology import LambdaCDM
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from astropy import constants as const
import astropy.units as u
from scipy.integrate import quad
import ezgal # BC03 model maker
import os
In [2]:
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7, Tcmb0=2.725)
In [3]:
# check to make sure we have defined the bpz filter path
if not os.getenv('EZGAL_FILTERS'):
os.environ['EZGAL_FILTERS'] = (f'{os.environ["HOME"]}/Projects/planckClusters/MOSAICpipe/bpz-1.99.3/FILTER/')
model = ezgal.model('bc03_ssp_z_0.02_salp.model')
model = model.make_exponential(1)
model.set_cosmology(Om=cosmo.Om0, Ol=cosmo.Ode0, h=cosmo.h, w=cosmo.w(0))
model.add_filter('g_MOSAICII.res', name='g')
model.add_filter('r_MOSAICII.res', name='r')
model.add_filter('i_MOSAICII.res', name='i')
model.add_filter('z_MOSAICII.res', name='z')
model.add_filter('K_KittPeak.res', name='K')
# Blanton 2003 Normalization
Mr_star = -20.44 + 5 * np.log10(cosmo.h) # abs mag.
# set the normalization
model.set_normalization('sloan_r', 0.1, Mr_star, vega=False)
In [4]:
# desired formation redshift
zf = 6.0
# fetch an array of redshifts out to given formation redshift
zs = model.get_zs(zf)
# Calculate some cosmological stuff
DM = cosmo.distmod(zs)
dlum = cosmo.luminosity_distance(zs)
Need to compute the cluster volume...
$M_{vir} = 4/3 \pi r^3_{vir} \rho_c(r<r_{vir}) = 4/3 \pi r^3_{vir} \Delta_c \rho_c$
if we let $\Delta_c = 200$ then
$M_{200} = 4/3 \pi r^3_{200} 200 \rho_c$ with $\rho_c = \frac{3H(z)^2}{8\pi G}$
or just $M_{200} = V_{200}200\rho_c$. So we'll make a function to calculate $\rho_c$. And we'll make use of the astropy units package to do all the unit analysis for us.
Don't forget that $H(z) = H_0E(z)$
The Schechter Function:
For Luminosity:
$\Phi(L) = \phi^\star \frac{L}{L_\star}^\alpha e^{-\frac{L}{L_\star}}$
For Magnitudes:
$\Phi(M) = \phi^\star\frac{2}{5}log(10) (10^{\frac{2}{5}(M_\star - M)})^{\alpha+1} e^{-10^{\frac{2}{5}(M_\star - M)}}$
In [5]:
def rho_crit(z, cosmo):
# convert G into better units:
G = const.G.to(u.km**2 * u.Mpc/(u.M_sun * u.s**2))
return 3 / (8 * np.pi * G) * cosmo.H0**2 * cosmo.efunc(z)**2 # Mpc^3
def schechterL(luminosity, phiStar, alpha, LStar):
"""Schechter luminosity function."""
LOverLStar = (luminosity/LStar)
return (phiStar/LStar) * LOverLStar**alpha * np.exp(- LOverLStar)
def schechterM(magnitude, phiStar, alpha, MStar):
"""Schechter luminosity function by magnitudes."""
# check to make sure things are all the same size
if isinstance(phiStar, float) and isinstance(magnitude, np.ndarray):
phiStar = np.ones_like(magnitude) * phiStar
if isinstance(alpha , float) and isinstance(magnitude, np.ndarray):
alpha = np.ones_like(magnitude) * alpha
if isinstance(MStar, float) and isinstance(magnitude, np.ndarray):
MStar = np.ones_like(magnitude) * MStar
MStarMinM = 0.4 * (MStar - magnitude)
return (0.4 * np.log(10) * phiStar * 10.0**(MStarMinM * (alpha + 1.)) * np.exp(-10.**MStarMinM))
In [6]:
def calc_Maglim(mlim, redshift, filter='i', cosmo=cosmo, model=model):
return mlim - cosmo.distmod(redshift).value - model.get_kcorrects(zf, filters=filter, zs=redshift)
def calc_phistar(redshift, cosmo=cosmo):
return 3.6 * cosmo.efunc(redshift)**2
def calc_alpha(redshift):
return -1.05 * (1 + redshift)**(-2 / 3)
def calc_fr(redshift):
return 0.8 * (1 + redshift)**(-1 / 2)
In [30]:
# So now we are going to calculate the volumes as a function of z
#M200 = mass_func(zarr) * u.solMass
M200 = 1e15 * u.solMass
V200 = M200/ (200 * rho_crit(zs, cosmo))
# Calculate the M_star values
Mstar = model.get_absolute_mags(zf, filters='i', zs=zs)
# calculate the abs mag of our limiting magnitude as a function of z
mlim = 23.5
#Mlim = Mstar - 2.5 * np.log10(0.4)
Mlim = mlim - cosmo.distmod(zs).value - model.get_kcorrects(zf, filters='i', zs=zs)
# Here are the Schechter function stuff from Liu et al.
phi_star = 3.6 * cosmo.efunc(zs)**2
alpha = -1.05 * (1 + zs)**(-2/3)
fr = 0.8*(1 + zs)**(-1/2)
In [8]:
LF = []
for phi, a, M_star, M_lim in zip(phi_star, alpha, Mstar, Mlim):
if M_lim < M_star - 2.5 * np.log10(0.4):
Mlimit = M_lim
else:
Mlimit = M_star - 2.5 * np.log10(0.4)
y, err = quad(schechterM, -30, Mlimit, args=(phi, a, M_star))
LF.append(y)
In [9]:
plt.figure()
plt.plot(zs, (LF * V200.value + 1) * fr)
ax = plt.gca()
ax.set_yticks(np.arange(0, 75, 10))
plt.xlim(0.1, 5)
plt.ylim(0, 80)
plt.xlabel('redshift')
plt.ylabel('N (r < r200)')
plt.grid()
In [10]:
calc_Maglim(23.5, 0.2)
Out[10]:
In [11]:
magarr = np.arange(-30, -14, 0.2)
In [12]:
redshifts = np.arange(0.2, 2, 0.2)
# make some colored lines for the plots
cm_subset = np.linspace(0.2, 0.8, len(redshifts))
colors = [cm.bone(x) for x in cm_subset]
The figure below is where we let both $\phi^\star$ and $\alpha$ evolve with redshift.
The vertical lines show the faint limit to our integral. In all cases we are integrating to $M=-30$ mag. So in the function above, we are integrating from the left to the right in this figure.
In [13]:
plt.figure(figsize=(8,6))
for z, c in zip(redshifts, colors):
plt.plot(magarr, schechterM(magarr, calc_phistar(z), calc_alpha(z), model.get_absolute_mags(zf, filters='i', zs=z)), c=c)
plt.axvline(calc_Maglim(23.5, z), c=c, label=f'{z:.2}')
plt.xlabel('Abs Magnitude')
plt.ylabel('$\Phi(M)$ [Ngal $Mpc^{-3}$]')
plt.semilogy()
plt.ylim(0.01, 20)
plt.legend(loc='upper left')
Out[13]:
The figure below is where we fix both $\phi^\star$ and $\alpha$ at some low redshift value, and do not allow them to evolve with redshift.
Here I've fixed them at $z=0.2$.
In my opinion, this looks a lot more like what we would expect for a LF. It still doesn't fully capture the evolution of the LF that I'd expect. I'd expect the
In [14]:
plt.figure(figsize=(8,6))
for z, c in zip(redshifts, colors):
plt.plot(magarr, schechterM(magarr, calc_phistar(0.2), calc_alpha(z), model.get_absolute_mags(zf, filters='i', zs=z)), c=c)
plt.axvline(calc_Maglim(23.5, z), c=c, label=f'{z:.2}')
plt.xlabel('Abs Magnitude')
plt.ylabel('$\Phi(M)$ [Ngal $Mpc^{-3}$]')
plt.semilogy()
plt.ylim(0.01, 10)
plt.legend(loc='upper left')
Out[14]:
In [52]:
LF = []
for phi, a, M_star, M_lim in zip(phi_star, alpha, Mstar, Mlim):
if M_lim < M_star - 2.5 * np.log10(0.4):
Mlimit = M_lim
else:
Mlimit = M_star - 2.5 * np.log10(0.4)
y, err = quad(schechterM, -30, Mlimit, args=(7, -0.8, M_star))
#print(M_star - M_lim, y)
LF.append(y)
Here is the above figure, with $\phi^\star$ and $\alpha$ fixed to their $z=0.2$ values.
In [53]:
plt.figure()
plt.plot(zs, (LF * V200.value + 1) * fr)
ax = plt.gca()
#ax.set_yticks(np.arange(0, 75, 10))
plt.xlim(0.1, 2)
#plt.ylim(0, 80)
plt.xlabel('redshift')
plt.ylabel('N (r < r200)')
plt.axhline(3)
plt.grid()
In [23]:
LF * V200.value + 1
Out[23]:
In [32]:
Mpiv = 6e14
zpiv = 0.46
fr2 = 0.68 * (M200.value / Mpiv)**-0.1 * ((1 + zs)/ (1 + zpiv))** -0.65
In [44]:
plt.plot(fr, fr2)
plt.plot(fr, fr3)
plt.plot([0,1], [0,1], c='k')
Out[44]:
In [41]:
from scipy.special import erf
N200 = LF * V200.value + 1
gz = 0.206 - 0.371 * zs
hz = -3.6 + 25.8 * zs
fr3 = gz * erf(np.log10(N200) - np.log10(hz)) + 0.69
In [42]:
fr3
Out[42]:
In [43]:
fr2
Out[43]:
In [54]:
calc_phistar(0.2)
Out[54]:
In [46]:
phi_star
Out[46]:
In [ ]: